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ABSTRACT 

We apply the MEGNO (Mean Exponential Growth of Nearby Orbits) technique to the dy- 
namics of Jovian irregular satellites. The MEGNO indicator is an effective numerical tool 
to distinguish between quasi-periodic and chaotic orbital time evolution for a given dynam- 
ical system. We demonstrate the efficiency of applying the MEGNO indicator to generate a 
mapping of relevant phase-space regions occupied by observed jovian irregular satellites. The 
construction of MEGNO maps of the Jovian phase-space region within its Hill-sphere is ad- 
dressed and the obtained results are compared with previous studies regarding the dynamical 
stability of irregular satellites. Since this is the first time the MEGNO technique is applied 
to study the dynamics of irregular satellites we provide a review of the MEGNO theory. We 
consider the elliptic restricted three-body problem in which Jupiter is orbited by a massless 
test satellite subject to solar gravitational perturbations. The equations of motion of the sys- 
tem are integrated numerically and the MEGNO indicator computed from the systems varia- 
tional equations. An unprecedented large set of initial conditions are studied to generate the 
MEGNO maps. The chaotic nature of initial conditions are demonstrated by studying a quasi- 
periodic orbit and a chaotic orbit. As a result we establish the existence of several high-order 
mean-motion resonances detected for retrograde orbits along with other interesting dynamical 
features. The computed MEGNO maps allows to qualitatively differentiate between chaotic 
and quasi-periodic regions of the irregular satellite phase- space given only a relatively short 
integration time. By comparing with previous published results we can establish a correlation 
between chaotic regions and corresponding regions of orbital instability. 

Key words: methods: n-body simulations — methods: numerical planets and satellites: ir- 
regular satellites — individual satellites: Carme - Ananke - Pasiphae - Themisto - Himalia - 
Sinope 



1 INTRODUCTION 

The existence of natural satellites in orbit around Solar System gi- 
ant planets has been known since Galileo discovered the four inner 
regular moons of Jupiter. Since then two classes of natural satel- 
lites have been identified to orbit the outer giant planets: regular 
and irregular satellites. The study of irregular satellites is interest- 
ing as they are located close to the Hill sphere boundary of the 
parent planet. Consequently Solar perturbations are present and 
chaotic behaviour is expected in the time evolution of their or- 
bits. Recent reviews on the subject of irregular satellites are given 
in Peale (1999); Jewitt & Haghighipour (2007); Nicholson et al. 
(2008). Regular satellites are characterised by circular orbits lo- 
cated close to the equatorial plane of the planet. This is the case 
of the Galilean satellites of Jupiter. The class of irregular satellites 
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differ in being on large eccentric and highly inclined orbits in either 
prograde (direct) or retrograde motion. 

The majority of the current population of giant planet irregu- 
lar satellites have been discovered by ground-based large aperture 
optical telescopes operating in dedicated survey programs over the 
past 10 years (Gladman et al. 1998, 2000, 2001; Sheppard & Jewitt 
2003; Holman et al. 2004; Sheppard et al. 2005, 2006). The most 
abundant irregular satellite population is observed to exist at Jupiter 
counting approximately 60 members at the current time. An inter- 
esting dynamical property of the phase space distribution of orbital 
elements is the clustering in distinct families for both prograde and 
retrograde irregular satellite members (Nesvorny et al. 2003, 2004). 

The regular Galilean satellites are suggested to have formed by 
accretion processes within a circumjovian planetary disk (Canup 
& Ward 2002; Mosqueira & Estrada 2003; Estrada & Mosqueira 
2006). Observed kinematic differences between regular and irreg- 
ular satellite orbits suggest a different formation mechanism of ir- 
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regular satellites. The most favoured formation scenario is capture 
from an initial heliocentric orbit. The key element of permanent 
capture is the necessity of a frictional force (due to a gaseous en- 
vironment) or energy transfer through close-encounter. Both pro- 
cesses are capable to dissipate orbital energy and provide a viable 
dynamical route to form a given population of irregular satellites. 
For a description of various capture scenarios we refer to Jewitt & 
Haghighipour (2007). 

Several authors have studied general stability properties of 
irregular satellites. Haghighipour & Jewitt (2008) studied the re- 
gion between Callisto and the innermost Jovian irregular satellite 
Themisto (30 Rj up < r < 80 Rj up where Rj up is Jupiters radius and 
r is the distance of the satellite to Jupiter). Despite observational 
evidence indicating that this region is void of a satellite population 
they showed that a large fraction of this region renders possible 
satellite orbits stable for at least 10 Myrs. Extensive stability sur- 
veys of irregular satellites close to the Hill sphere of the giant plan- 
ets have been carried out by Carruba et al. (2002); Nesvorny et al. 
(2003) by recording the lifetimes of irregular satellite test particles 
in various parameter surveys. Yokoyama et al. (2003) studied the 
region 250 Rj up < r < 370 Rj up for prograde Jovian irregular or- 
bits in the semi-major axis and inclination plane finding evidence 
of the presence of secular perturbations (Whipple & Shelus 1993; 
Cuk & Burns 2004; Nesvorny & Beauge 2007). Stability properties 
of irregular orbits at the outer regions around the giant planet's Hill 
sphere and beyond were studied by Shen & Tremaine (2008). 

Motivated to study the phase space topology structure of ir- 
regular satellites in detail we applied the MEGNO chaos indi- 
cator (Cincotta & Simo 2000; Gozdziewski 2001; Gozdziewski 
et al. 2002) to qualitatively differentiate between quasi-periodic and 
chaotic phase space regions. Initial tests showed that MEGNO is 
efficient in showing chaotic regions using relative short integration 
times of the orbit. In this work we outline the basic principles of the 
MEGNO indicator as this is the first time this technique is applied 
to the dynamics of irregular satellites. 

The structure of the paper is as follows. In section 2 we present 
the model and numerical methods used as well as initial conditions 
and definitions of angular variables. Section 3 presents a review of 
the MEGNO chaos indicator and an outline of its computation. Its 
relation to the Lyapunov exponent is reviewed and tests on numer- 
ical accuracy are presented. Section 4 describes the construction 
of dynamical maps of Jovian irregular satellites as studied in this 
work. Section 5 presents and outlines essential steps to obtain the 
secular system of the time variation of a satellite's Keplerian ele- 
ments using a time-averaged running window. Section 6 presents 
our results with comparison to previous work and section 7 con- 
cludes this work. 



2 MODEL, NUMERICAL METHODS AND INITIAL 
CONDITIONS 

The results obtained in this work are based on the elliptic restricted 
three-body problem. We integrate the equations of motion of the 
system (Morbidelli 2002) 

d 2 rj _ k 2 m A 2 | Tj-Tj rj \ 

dt 2 ~ |nP n+ .LJ^^rj-rtf \rj\*)' W 

Here n = 2 in a jovian centric reference frame with mo denoting 
the mass of Jupiter and k 2 denotes the Gauss gravitational constant. 
The positions (relative to Jupiter) and masses of the satellite and 
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Figure 1. Graphical illustration of the relationship between angles in a he- 
liocentric reference system and a planetocentric system. In both panels T 
denotes the reference direction from which inertial angles are measured and 
the orbital motion is counterclockwise. Panel A: Shows the orbit of (J)upiter 
with respect to the (S)un. Jupiters orbit is strongly exaggerated in its eccen- 
tricity for reasons of clarifications. In the current situation the longitude of 
pericenter of Jupiter is toJ e/ = 135°. Panel B: The orbit of (S)un with respect 
to (J)upiter. Here the longitude of perijove U5 J ° V = 315°. 



the Sun are (r\,m\) and (^rni), respectively. In this model the 
satellite is orbiting Jupiter and its orbit is subject to Solar pertur- 
bations. In all integrations we regard the satellite as a test particle 
of zero mass. Oblateness effects from Jupiter are not considered in 
this model and perturbations from other planets are omitted as well. 
We consider this simplified model for two reasons. First the pur- 
pose of this paper is to demonstrate and apply the MEGNO chaos 
indicator to the dynamics of irregular satellites and secondly we 
aim to qualitatively identify chaotic regions originating from Solar 
perturbations only. If we were to include oblatenes and planetary 
perturbations cause and effect would be very difficult to isolate. 
The orbit of the satellite was obtained by the numerical integra- 
tion of Eq. (1) using the 1) 15th-order Radau integration algorithm 
as implemented in the latest version of the MERCURY package 
(Chambers & Migliorini 1997; Chambers 1999) and 2) the Gragg- 
Bulirsch-Stoer (GBS) extrapolation algorithm (Hairer et al. 1993) 
as implemented in the CS -MEGNO code (Gozdziewski 2001). 

The Radau integrations were used to study individual satel- 
lite orbits on the long-time scales and the GBS integrator was used 
for the generation of MEGNO maps over the observed irregular 
satellite phase space region. Initial conditions (geometric cartesian 
elements relative to Jupiter's centre of mass) for the Sun and ob- 
served irregular satellites have been obtained from the JPL Horizon 
1 Ephemeris generator (Giorgini et al. 1996) at the epoch 01-Jan- 
2005 (12:00 UT). It is worth to note that the Horizon Ephemeris 
generator cannot output osculating Keplerian elements of the Sun 
with respect to another object in the Solar System. Only carte- 
sian vectors can be retrieved for the Sun relative to Jupiter. For the 
computations of the stability maps (see section 6.1) we have trans- 
formed the Sun's cartesian elements to Keplerian elements using 
the MCO_x2el.F subroutine as implemented in the latest version of 
MERCURY. The transformation introduces round-off errors not 
larger than 10 -10 in the absolute errors which is much smaller than 
the error of observed elements of irregular satellites (e.g 15 meters 
in the semi-major axis). The retrieved initial conditions are referred 
to the ecliptic and mean equinox (xy-plane is the orbit of Earth) at 
reference frame ICRF/J2000.0. In this work we will present our re- 
sults in a planetocentric reference system where Jupiter is at the 
centre and the ecliptic is the reference plane. Thus we denote the 
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planetocentric elements 2 of a satellite as (a, e,I, co, £l,M) (where M 
is the mean anomaly) and we use the subscript to indicate the 
Sun's orbit relative to Jupiter and the subscript / to denote the orbit 
of Jupiter in the heliocentric reference system. For the longitude of 
perijove of prograde irregular satellites we use the usual definition 
05 = £1 + CO. For retrograde satellites we use 05 = £1 — CO. The mean 
longitude of an irregular satellite is then given by X = 05 + M. 

In order to compare our data with previous published results 
in the literature using both planetocentric and heliocentric orbital 
elements we briefly point out the relationship of angles measured 
in the two reference frames. Changing the coordinate system from 
heliocentric to a jovicentric system leaves the semi-major axis and 
eccentricity unchanged as these quantities are invariable under a 
coordinate transformation. In Fig. 1 we demonstrate the relation- 
ship between the longitude of pericenter of the two bodies in the 
two reference system. If USj el is the longitude of Jupiter in the 
heliocentric system and G5q v denotes the longitude of the Sun in 
the jovicentric system, then from geometric arguments we have 
^hel _ Q^jw + jgQo ^ s j m ji ar ar g Umen t w iH lead to the relation- 
ship X h / 1 = X J ° V + 180° relating the mean longitudes of Jupiter and 
the Sun. The orbital inclination (/©,//) and the argument of node 
(£Iq,£Ij) are unchanged under the transformation. 



3 THE MEGNO CHAOS INDICATOR 

The time evolution of irregular satellite orbits exhibits both quasi- 
periodic and chaotic dynamics (Whipple & Shelus 1993; Saha & 
Tremaine 1993; Goldreich & Rappaport 2003). An efficient nu- 
merical method to detect phase space regions resulting in chaotic or 
quasi-periodic initial conditions is provided by the MEGNO (Mean 
Exponential Growth of Nearby Orbits) factor (or indicator). The 
MEGNO technique was first introduced by (Cincotta & Simo 2000; 
Cincotta et al. 2003) and were originally inspired from the concept 
of 'conditional entropy of nearby orbits' (CENO) (Cincotta & Simo 
1999; Gurzadyan 2000). It can be applied to any dynamical system 
with more than two degrees of freedom and has found widespread 
applications in dynamical astronomy ranging from galactic dynam- 
ics to stability analysis of extrasolar planetary systems and Solar 
System small body dynamics (Cincotta et al. 2003; Gozdziewski 
2001, 2002, 2004; Breiter et al. 2005). 

By numerically evaluating the MEGNO factor after a given 
integration time one obtains a quantitative measure of the degree of 
stochasticity of the system. One method (Morbidelli 2002; Dvo- 
rak et al. 2005, and references therein) to discriminate between 
ordered (or regular) and chaotic (very often associated with or- 
bital instabilities) satellite orbits, is the calculation of the system's 
Lyapunov characteristic exponents (LCE) or Lyapunov character- 
istic numbers (LCN). From the theory of dynamical systems each 
system has a spectrum of Lyapunov exponents (possibly complex 
eigenvalues) each associated to a given eigenvector of the system. 
The Lyapunov exponents describe the rate of change of its cor- 
responding eigenvector in time. The number of positive (or van- 
ishing) Lyapunov exponents indicates the number of independent 
directions in phase space along which the satellite orbit exhibits 
chaotic (or quasi-periodic) behavior (Morbidelli 2002). In partic- 
ular, MEGNO is closely related to the maximum Lyapunov expo- 
nent (MLE or sometimes referred to maximum Lyapunov numbers, 



MLN) providing a quantitative measure of the exponential diver- 
gence of nearby orbits and belongs to the class of fast Lyapunov 
indicators (Morbidelli 2002). 

In this work, we apply the MEGNO criterion to the dynamics 
of Jovian irregular satellites. Details on the MEGNO concept and 
its numerical computation can be found in Cincotta et al. (2003) 
and Gozdziewski (2001). Since the MEGNO technique is applied 
for the first time to the dynamics of irregular satellites, we give a 
short review of the most important aspects of MEGNO and provide 
additional information on its numerical computation. 

3.1 The Lyapunov Exponent 

The maximum Lyapunov exponent y, provides a useful quantitative 
measure to study the dynamical nature of the time evolution of a 
satellite orbit in phase space and is defined as 

t^t-tQ \ 8(/ ) / t^°°t-t Jo 8(s) 

where 5(?) is the variational vector and measures the distance in 
phase space between two initially nearby orbits as a function of 
time. The second equality is easily realised by change of variables 3 . 
For y > 0, an initial separation grows exponentially in time (defi- 
nition of chaotic motion) at the rate or decays, if y < 0. In the 
case y = 0, the time rate of change of the variational vector 5 is 
0. Since the (elliptic) restricted three-body problem is a conserva- 
tive system the case y < is never encountered and would indicate 
the presence of dissipative forces. In practical computations only 
a finite time estimate y(t) of y(s) is obtained after integration time 
t. In our computations, we will pay some attention on the proper 
choice of t. Usually the convergence of y(t) is slow and a reliable 
numerical estimate of y requires a long integration time of the dy- 
namical system. When exploring the dynamics of a large portion 
of phase space, short integration times are preferred when explor- 
ing the phase-space structure. For classic computations of the Lya- 
punov exponents from the variational equations we refer to (Benet- 
tin et al. 1980; Wolf et al. 1985). 

3.2 Mathematical properties of MEGNO 

A faster convergence property is obtained from the mean exponen- 
tial growth factor of nearby orbits. The MEGNO indicator is closely 
related to the definition of y and is defined as 

Y ® = lfw{ sds > (3) 
t Jo o{s) 

along with its time-averaged mean value 

(Y)(t) = \ fY{s)ds. (4) 
t Jo 

Cincotta & Simo (2000) showed that (Y)(t) converges faster to its 
limit value than the Lyapunov characteristic number. In the former 
definition the relative rate of change of the separation vector 5/5, 
is weighted with time during the integration giving preference to 
the memory of late evolutionary behavior of the separation vector 
(Morbidelli 2002). The time- weighting factor introduces an ampli- 
fication of any stochastic behavior, allowing an early detection of 
chaotic motion (Gozdziewski 2001). 



2 the orbital inclination of a satellite is measured relative to the ecliptic. 



3 If y = 8(*), then f>(d(s)/5(s))ds = fl(l/y)dy = \n(y(t)/y(0)) 
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Figure 2. Results from 2-body Kepler integrations. Right panel: A) (RTOL,ATOL) = (5.0 • 10" 11 , -lO" 10 ), B) (RTOL,ATOL) = (5.0 • lO" 13 , -lO" 10 ), 
C) (RTOL,ATOL) = (1.0 • KT 11 , -lO" 10 ), D) (RTOL,ATOL) = (1.0 • 10" 10 , -lO" 10 ), E) (RTOL,ATOL) = (5.0 • KT 11 , -lO" 10 ), F) (RTOL,ATOL) = (1.0- 
10- 16 ,1.0-10- 16 ),G) (RTOL,ATOL) = (1.0- 10~ 15 , 1.0- 10~ 15 ), H) (RTOL,ATOL) = (1.0 • 10" 14 , 1.0 • 10" 10 ), I) (RTOL,ATOL) = (5.0 • 10" 15 , 1.0 • 10" 12 ). 
The area marked with vertical lines corresponds to case D and horizontal line corresponds to case B. The final value after 10 million years of (Y) for case F) 
was 1 .999757, case G) 2.000367, case H) 1 .636787 and case I) 1 .798234. 



Computing the time evolution of Eq. (4) for a set of initial 
conditions allows the study of the dynamical properties in a given 
phase space region. Following Cincotta & Simo (2000); Cincotta 
et al. (2003), if the motion is of quasi-periodic nature, then 5 grows 
linearly in time and (Y) will exhibit asymptotic oscillations about 
2 for t — > oo. In the case of chaotic initial conditions, (Y)(t) =yt/2 
for t — > oo. The latter equation relates the Lyapunov characteristic 
number to the limiting value of (Y) (t) at integration time t. A linear 
best fit to (Y)(t) would recover the Lyapunov characteristic num- 
ber (Gozdziewski 2001) . In summary, the time-averaged MEGNO 
((Y)(t)) converges faster to its limit value compared to the stan- 
dard calculation of the Lyapunov characteristic number. This prop- 
erty allows a more rapid exploration of the dynamical phase space 
structure. 



the first term describes the variation in the 2-body Kepler orbit and 
the variation in perturbations (or interactions) is 
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Here, = rj — r/, Sr i; - = 8r ; — 8r ; and broj = 8r 7 — Sro, where 
the zero superscript denotes the central body. Eqs. (6)-(7) are com- 
puted in a straightforward way once the initial conditions have been 
defined for the variationals and the initial osculating orbits. 



3.4 Numerical accuracy tests 



3.3 Numerical calculation of Y(t) and (Y) (t) 

In practice, Y(t) and (Y)(t) are calculated by rewriting Eqs. (3) and 
(4) into two differential equations (Gozdziewski 2001) 



dx 8 , dw „ x 

— = ^t and — = 2 -, 
dt 8 dt t 



(5) 



where Y(t) = 2x(t)/t and (Y)(t) = w(t)/t. To obtain x(t),w(t) 
the two first-order coupled differential equations are solved along 
with the equations of motion. Initial condition for the variationals 
are chosen using a random generator. The quantities 5 = / 5r and 
5 = / 5v of the /'th body are obtained by solving the variational 
equations 8 = bv and 8 = bv (Mikkola & Innanen 1999). At each 
time step we have for the i— th body the variationals 



5r ; 



5v; 



(6) 



5, = -* {m + m) [**!^\-^ (7 ) 



The MEGNO indicator is numerically computed as outlined in the 
previous section using the Gragg-Bulirsch-Stoer integration algo- 
rithm (Hairer et al. 1993). All computations are carried out using 
double precision arithmetic. To control the numerical errors during 
computations two accuracy parameters are necessary for a GBS 
integration. The two parameters control the absolute (5) and rela- 
tive (e) error tolerances for any given integration and both needs 
to be specified for a given accuracy requirement. In practical com- 
putations the usual choice of (8,8) falls into the range from one 
part in 10 9 down to the limit of the machine precision (10 -16 ) 
(depending on architecture). To find a suitable (8, e) set, we inte- 
grated several 2-body Kepler problems for 10 million years by fol- 
lowing the orbit of Himalia. In that case the motion is expected to 
be quasi-periodic, and we use it as a test case to determine confi- 
dence of the computed MEGNO indicator which should converge 
to (Y)(t) = 2.0. Initial numerical tests showed that the convergence 
of (Y) depends on the accuracy of the numerical integration. This 
property was already outlined in Gozdziewski (2001). We inte- 
grated Himalia' s orbit (without Solar perturbations) for different 
combinations of the tolerance parameters. In total we considered 
(14 x 14) combinations with (8,8) ranging from 10~ 10 to 10 -16 . 
The result of the integrations is shown in Fig. 2. The left figure 
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Figure 3. Time evolution of Y(t) (solid curve) and (Y)(t) (dotted curve) 
from a numerical integration over 1 Myrs for the retrograde satellite Ananke 
including Solar perturbations. The solid horizontal line indicates (Y) = 2.0. 
The upper time line expresses integration time in units of Jupiters orbital 
period Pj w 12.5 yrs. 



panel shows the numerical value of (7) (gray- scale color coded) 
at the end of the integration for a given combination of (5,8). The 
right panel in Fig. 2 shows the time evolution of (Y)(t) of a few 
selected combinations of the tolerance parameters. The (5, e) com- 
binations resulting in the cases denoted by (A-E) are evidently of 
poor numerical accuracy showing (Y) ^> 2.0 after 10 million years. 
We interpret this result as artificial numerical chaos, the result of 
the poor resolution of time descretisation in the integration algo- 
rithm. The cases (F,G,H,I) show a significant better agreement with 
the expected value of (7) = 2.0 with some numerical fluctuation 
of (Y) observed during the integration. From the numerical experi- 
ments we fix (5,8) = (10" 16 , 10 -15 ) for all MEGNO computations 
in the present work. In addition, we monitor the relative energy er- 
ror dE/E and note the maximum value reached during a given in- 
tegration. For our choice of the tolerance parameters the average 
maximum relative energy error is on the order of dE/E ~ 10 -11 
over 1 million years. By examining several spot tests on the time 
evolution of the relative energy error no systematic trends were ob- 
served only exhibiting a random walk over time. For all orbits the 
maximum relative energy error were smaller than 10 -12 though in 
the restricted three-body problem this does not reflect the accuracy 
of the orbit of a massless test satellite. For that reason we follow the 
same approach as outlined in Nesvorny et al. (2003) and monitored 
the Jacobi constant during the numerical integration and found a 
maximum relative change of this quantity to be on the order of typ- 
ically 10 -6 . We also computed the absolute error of the semi-major 
axis (|fl©(f) — a© |) of the Sun and found it to be less than a few 
parts in 10 -12 . This means a preservation of 8 significant digits of 
the semi-major axis over 60000 years. Since we are not aiming at 
generating high-accuracy ephemerides of irregular satellite orbits 
we consider these tests as sufficiently reliable in order to establish 
confidence in our results presented in this work. In addition, un- 
certainties in observed orbital elements of irregular satellites are 
assumed to be much greater than rounding/truncation errors intro- 
duced by the numerical integration algorithm. 




Figure 4. Time evolution of (Y) for different selected irregular satellites 
over a time scale of 1 million years. The shown satellites are Th(emisto), 
Hi(malia), Ca(rme) and Si(nope). 

3.5 Numerical examples of test orbits and limitations 

As an example Fig. 3 shows the time evolution of Y(t),(Y)(t) 
for the irregular satellite Ananke over a time span of 1 Myrs. In 
this integration, gravitational perturbations from the Sun were in- 
cluded and it is observed that Ananke's orbit exhibits a weak sign 
of chaoticity. It is evident that no clear convergence of (7) ^2 
is observed. Also the time evolution of Y(t), (Y)(t) suggests that 
the time-averaged MEGNO ((F)) is more reliable to indicate the 
presence of chaotic behavior as Y(t) exhibits large variations about 
Y(t) =2.0 during the integration. 

In Fig. 4 we show the results of calculating (7) of four se- 
lected irregular satellites (Carme, Himalia, Sinope and Themisto) 
over 1 Myr. Our results suggest that both Themisto and Himalia 
(both on prograde orbits) shows chaotic behaviour as their calcu- 
lated (Y) show clear sign of divergence from 2.0. The contrary is 
seen for the orbits of Carme and Sinope (both on retrograde orbits) 
indicating convergence (although slow for Carme) towards 2.0 in- 
dicating quasi-periodic time evolution. This result is in agreement 
with what is expected from numerical simulations. Prograde orbits 
tend to be less stable compared to retrograde orbits (Hamilton & 
Krivov 1997; Nesvorny et al. 2003). It is noteworthy that although 
the calculation of (7) is a powerful numerical tool for the detec- 
tion of chaotic initial conditions it is necessary to be cautious when 
interpreting results. As mentioned previously two of the irregular 
satellites (Themisto and Himalia) show a weak sign of chaotic be- 
haviour with (7) diverging from 2.0 with | (7) - 2.0| < 0.0003 after 
200000 years. In a MEGNO map such an initial condition would 
appear as quasi-periodic considered the large range in (7) cover- 
ing several orders of magnitude. Therefore it is imperative to point 
out that every numerical tool used to differentiate between quasi- 
periodic and chaotic behaviour is only capable of showing quasi- 
periodicity up to the integration time. On the contrary once chaotic 
behaviour (excluding numerical chaos) has been detected the or- 
bit can be claimed chaotic for all times. If each initial condition in 
the MEGNO maps presented in this work were calculated on time 
scales similar to the age of the Solar System it is likely that no 
quasi-periodic orbits are detected. In our maps we can with high 
confidence associate chaotic regions with unstable orbits where the 
satellite either escapes or experience a collision with one of the 
inner moons or if the eccentricity grows high enough (due to the 
Kozai mechanism) the satellite may even collide with Jupiter itself. 
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Figure 5. Demonstrating the process of successive smoothing using a running window time average applied on the semi-major axis. See text for details. 



4 CONSTRUCTION OF MEGNO, ME, MS AND MI MAPS 

In this work we color code the MEGNO indicator on a 2- 
dimensional phase space section mapping out the dynamics 
of irregular satellites. In particular we show the dynamics in 
(a J)— space, where a is the semi-major axis and / denotes the orbit 
inclination. We chose initial conditions in which orbit inclinations 
are relative to ecliptic (see IC section). For a given range in semi- 
major axis and inclination the grid of initial conditions in the maps 
are given by 

. Aa . Umax ~ a min . . irv . 

ai = a min + — i = a min H i (10) 

N x N x 

and 

T T AI . I max ~ hnin . 

Ij = Imin + — J = Imin H 7, (1 1) 

y y 
where i = 0, ...N x and j = 0, ...N y are integers defining the 
grid resolution. Chosing a large (N x ,N y ) will result in a more 
detailed mapping of phase-space structures in (a, I) space. The 
range in semi-major axis is chosen to span a e [0.04,0.20] 
AU ([0.11,0.56] R H , [83,419]R Jup ). Satellite inclinations cover the 
range / e [0, 180] degrees considering both prograde and retrograde 
satellite orbits. In all MEGNO maps we chose to color code initial 
conditions resulting in quasi-periodic motion by blue (see the elec- 
tronic version of this work). In the quasi-periodic case we have 
— 2.0 1 < 0.01 at the end of the numerical integration. In ad- 
dition, we also provide maps showing the maximum eccentricity 
(ME) of an orbit at a given initial grid point. ME maps were con- 
structed in parallel with the MEGNO calculations. The osculating 
elements were measured at each completed Jovian period and the 



maximum value is determined by comparing with a previous mea- 
surement of the eccentricity. 



5 SMOOTHING ORBITAL ELEMENTS 

To study the secular system of irregular satellites we have to filter 
out the fast frequencies. Saha & Tremaine (1993) already pointed 
out that the unfiltered variations in the action elements (semi-major 
axis, eccentricity and inclination) are larger than the filtered time 
series of a given element. This means that the fast variations (high- 
frequency terms) are much larger in amplitude than the slow vari- 
ations (long-period terms). Hence the secular system is masked by 
the high frequencies. In this work an initial preliminary study of 
several test orbits confirmed this dynamical behaviour and we find 
that the most interesting dynamical features are to be found in the 
secular system. To obtain the secular system one can either average 
out the fast frequencies by applying a digital filter in either the time 
or frequency domain (Carpino et al. (1987); Quinn et al. (1991); 
Saha & Tremaine (1993); Michtchenko & Ferraz-Mello (1995)) or 
by simply averaging out all quasi-periodic oscillations with a run- 
ning window average (Morbidelli 1997; Morbidelli & Nesvorny 
1999). 

In order to study the secular system we smooth the orbital 
elements of a given time sequence Aj using the SMOOTH function 
as implemented in IDL 4 . The smoothing procedure is a running 

4 IDL stands for Interactive Data Language. For more information 

http : //www. ittvis . com/Product Services /IDL . aspx 
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Figure 6. Time evolution of orbital elements for initial condition IC-I (left panel, chaotic) and initial condition IC-II (right panel, quasi-periodic). The full 
data set (panel a, c, e and g) and corresponding smoothed data set (panel b,d,f and h) are shown in each panel. From top to bottom panel: Semi-major axis 
(a+b). Eccentricity (c+d). Inclination (e+f) and the angle 03 - G3 (g+h). The smoothing window applied to the data set in both panels are identical. Only the 
initial conditions has changed. Half of the window widths for a given element can be seen as a black bar at the beginning and end of the integration time. 
Each initial condition has been integrated for 60000 years corresponding to the integration length in the MEGNO maps. The window width are w/2 x AT = 
2500/2 x 40 days = 137 years for panels b,d,g and w/2 x AT = 20000/2 x AT = 1095 years for panel f. Note that 03© = 05/ - 180°. 



window average applied on the full data set as obtained from a 
numerical integration. The output of a given numerical integration 
is given as a time sequence of orbital elements sampled at regular 
intervals of length AT. For a given sequence of an orbital element 
A( = A(ti) with ti = t$ + iAT (for i = 0, 1 , . . . , N) the smoothed 
(secular) sequence R( is given by 

( 1 yW-\ a . • _ W-l at W+l 

Ri=\ ™H=Q A i+i-{*?Y i--2",...,JV -J-. 
1 Aj otherwise. 

Here N is the number of data points from the numerical integration 
and w is the (running) window width over within which the original 
data set is averaged on. In units of time the window width is simply 
wAT. The window width represents a free parameter and basically 
controls the supression of dynamical features seen in the data set: a 
too short window will have little smoothing effect and thus the fast 
frequencies are retained; a too large window width will suppress 
long-period features that appear on secular time scales. In practice 
some experimentation is needed in order to determine a satisfying 
window width. In this study we have performed an extensive survey 
of various window widths. In each plot showing smoothed orbital 
elements we have chosen the most appropriate window width in 
order to highlight the secular changes of a particular orbit. 

Some care has to be taken when smoothing angular vari- 
ables. As already mentioned by Saha & Tremaine (1993) filtering 
(or smoothing) an angular variable is problematic as angles may 
change discontinuously over the time sequence (i.e from — n to it). 
Applying a running window smoothing procedure directly to an 



angular variable introduces spurious results around discontinuities. 
To circumvent this problem, we transform a given angular variable 
to a continuous signal. Let 9 be an angular variable changing dis- 
continuously at times during the numerical integration. Then we 
transform to the following continuous variables 

p = Acos(9) (12) 
q = -Asin(e), (13) 

where A is some suitable constant. Here we chose A = 1. The 
IDL SMOOTH smoothing procedure is then applied to each of 
the quantities p,q and we obtain p,q after which we obtain the 
smoothed angle 9 from 9 = arctan(g/ p). 

An important note about the SMOOTH routine is the follow- 
ing. The averaging behaviour at the beginning and end of the orig- 
inal time sequence depends on the optional edge .truncate key- 
word passed to the SMOOTH function. If this keyword is enabled 
the smoothing procedure might introduce false/misleading aver- 
ages at the beginning and end of the smoothed signal. Details can be 
found in the IDL documentation. In addition if this keyword is dis- 
abled then it is important to note that Ri = A[ for the first data points 
up to (but not including) (w — 1 ) /2 and R[ = A/ from N—(w+l)/2 
(but not including) to N. 

We demonstrate the effect of successively applying time- 
averaging smoothing windows to the time evolution of the semi- 
major axis in Fig. 5. Fig. 5a shows the whole signal over 60000 
yrs with a sampling frequency of 40 days. At this stage a secu- 
lar period of about 2400 years is clearly present in the frequency 
spectrum of the signal. Fig. 5b shows the first 200 years of the full 
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Figure 7. MEGNO maps in (a,/) -space of test satellites with e = 0.20, (Do = 0°,£2 = 0°,M = 90°. The MEGNO indicator (Y) is color coded on a linear 
scale from 1.5 to 4 with (Y) = 2 indicating quasi-periodicity and (Y) = 4 denoting chaotic initial conditions. See the electronic version for colors. In both panels 
the symbols denote locations of irregular satellites. Themisto is at 0.05 AU (0.14 Rh, 105 Rj up and indicated by a x symbol). Himalia family = (+)-symbols 
and Carpo = (triangle) at 0.112 AU (0.32 R H ,235 Rj up ). We plot family members for which \e(t) - 0.20| < 0.05. Ananke (black • symbols), Carme (pink 
rhomb symbols) and Pasiphae family (magenta square symbols). Left panel a: Three regions showing interesting dynamical features are shown by rectangles. 
Right panel b: High resolution zoom of the upper right region shown in panel a. Arrows indicate the locations of retrograde mean-motion resonances. Nominal 
locations in are given in Table 1. Also shown are the locations (IC-I and IC-II) of two test orbits. As already detected by Yokoyama et al. (2003) note the small 
stable region at (a Jo) = (0.142 AU, 10°). 



data set (thin line) which is dominated by the orbital frequency of 
the satellite and an approximately 12 year period (Jupiter's orbital 
period). When applying two successive smoothing windows to the 
raw data with window width 3 yrs and 12 yrs we obtain the aver- 
age signal overplotted as a thick line in Fig. 5b. Fig. 5c shows a 
zoom of the smoothed signal over the whole integration time (note 
the difference in range of the semi-major axis). Fig. 5d (thin line) 
shows the first 3500 yrs of the previous signal. This time we detect 
a 34 year periodicity in the semi-major axis. Applying a smooth- 
ing window removes the 34 yr period (thick line in Fig. 5d) and 
the smoothed signal over the 60000 yrs is shown in Fig. 5e. Fur- 
thermore a 140 yr period signal is present in the semi-major axis as 
shown in Fig. 5. Applying a fourth smoothing window to the orig- 
inal raw signal now also removes this period and we end up with 
the secular signal (2400 yr period) shown in Fig. 5f. 

In Fig. 6 we show examples of the effect of the smoothing 
procedure when applied to the numerical solution of two different 
initial conditions located close to each other in (a,/)-space. In both 
panels we show the semi-major axis, eccentricity, inclination and 
the angle 05 — G5 for each initial condition. The left (right) panel is 
the result of integrating initial condition IC-I (IC-II) as indicated in 
Fig. 7. Each orbit were integrated for 60000 years, the integration 
length of each grid point in the MEGNO maps. The black bars (or 
sometimes scattered data points) shows half the window width (see 
figure caption for more details) and has been experimentally deter- 
mined based on qualitative judgment with the goal to enhance the 
underlying dynamical changes in a given osculating element. As 
already pointed out in Saha & Tremaine (1993) we again observe 
that most of the variation in the orbital elements are found in the 
fast frequencies. A visual inspection of the time evolution of initial 
condition IC-I confirms its chaotic nature as correctly identified by 
MEGNO. For comparison a time-running window with identical 



window width has been applied to initial condition IC-II. Its quasi- 
periodic nature is clearly visible over the considered time span. 
From experimentation with the window width we observed that in- 
creasing the width averages out more and more quasi-periodic os- 
cillations. No chaotic behaviour has been observed when increas- 
ing the window width. Following (Morbidelli & Nesvorny 1999, 
p.301) by applying a time-running window smoothing procedure to 
a given osculating element we obtain the corresponding proper el- 
ement. If the orbit evolves quasi-periodically (regular motion with 
a limited number of frequencies) in time then the corresponding 
proper element is a constant of motion. On the other hand if the 
corresponding proper element is chaotic then it exhibits a random 
walk in time. The chaotic nature of integrating initial condition IC- 
I particularly manifests itself in the time evolution of 05 — G5 as 
shown in Fig. 6h (left panel). This angle repeatedsly changes from 
libration to circulation; this is characteristic of chaotic behaviour 
(i.e perturbed pendulum model). For comparison this angle circu- 
lates for initial condition IC-II (quasi-periodic) over the entire inte- 
gration time span without any change in oscillation mode. 

6 RESULTS AND DISCUSSION 
6.1 MEGNO, a max and e max maps 

Our results on computing the MEGNO indicator over a large grid 
in (ao,A))-space is shown in Fig. 7. The figure shows the grid 
(aoM e ([0.04,0.20] AU, [0.0°, 180°] region. The observed pop- 
ulation of prograde and retrograde satellites are shown by vari- 
ous symbols (see figure caption for details). We obtained the os- 
culating elements of the irregular satellites from the JPL Hori- 
zon Ephemeris system (see earlier section on initial conditions). 
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Figure 8. Maximum semi-major axis (a max ) of a given initial condition in (a, I) -space for a irregular test satellite. Left panel a: Same range in semi-major axis 
and eccentricity as for the MEGNO map with the maximum semi-major axis color coded in the range [0.04,0.20] AU ([0.11,0.56] Rh, [84,419] Rj up ). Right 
panel b: Zoom of the rectangular area as indicated in the left a-panel figure. 



We chose to consider this part of phase-space in order to compare 
with previous results as published by Carruba et al. (2002) 5 and 
Yokoyama et al. (2003). These authors explore similar phase space 
regions although they consider a much smaller grid of initial con- 
ditions. Reproducing previously published results motivated us to 
apply and compute the MEGNO indicator over a much larger re- 
gion occupied by the outer irregular Jovian satellites. In our work 
at initial time the eccentricity was set to eo = 0.20 and the remain- 
ing Kepler elements (coo,^o) were set to zero with Mq = 90°. 
We calculated the MEGNO indicator for 35100 initial conditions 
(N x ,N y ) = (195, 180). The orbit of each initial condition were inte- 
grated for 5 x 10 3 Pj years (Pj is the orbital period of Jupiter). 

To save computing time we stopped a given integration as 
soon as (Y) ^ 10.0. The MEGNO indicator is then color-coded 
with (7) =2.0 indicating quasi-periodic motion and we only plot 
(Y) ^ 4 to enhance the contrast of dynamical features in the tran- 
sition region were the dynamics change from quasi-periodic to 
chaotic. Figure 7a shows several interesting regions of chaotic 
and quasi-periodic nature with unprecedented detail. We decided 
to compute high resolution maps to study this interesting region. A 
zoom plot is shown in Fig. 7b and corresponds to the area within 
the box in the upper right corner of Fig. 7a. The labels IC-I and IC- 
II correspond to initial conditions for two hypothetical (retrograde) 
irregular satellites. 

In parallel with computing (Y) to detect chaos we also 
recorded the maximum values of the test satellites semi-major axis 
and eccentricity every orbital period of the Sun. Fig. 8 and 9 show 
the corresponding maps. It is now apparent that the global chaotic 
region results in either escape or collisions after only 5 x 10 3 Pj 
years. From this result we can conclude that the large chaotic re- 
gion detected in the MEGNO map is strongly correlated with ei- 
ther escape from the Jovian Hill sphere or collisions with Jupiter 

5 The authors uses different units for the semi-major axis. Carruba 
et al. (2002) measures distance in units of Jupiter's Hill radius Rh, and 
Yokoyama et al. (2003) measures distance in units of Jupiter's radii Rj up . 
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Table 1. Nominal locations of retrograde MMRs as computed from 
Eq. (14). The masses of Jupiter and the Sun and the semi-major axis of 
Jupiter are taken from Murray and Dermott (Table A.l and A.2). 



itself. It is also interesting to note that over a large range the orbit 
size (semi-major axis) remains close to its initial value. Further- 
more we also note that the locations of mean motion resonances 
are not visible in the zoom plot (Fig. 8b) as otherwise indicated 
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Figure 9. Same as Fig. 8 but using e max for the colorcoding. Note the horizontal feature of enhanced eccentricity for both pro- and retrograde satellites. This 
feature was also reported to exists from a 1 Myr integration (Yokoyama et al. 2003, fig.6). 



by MEGNO. This indicates that the semi-major axis is unaffected 
by the presence of mean motion resonances and MEGNO is effi- 
cient in detecting mean motion resonances. A different picture is 
obtained when studying the maximum eccentricity attained during 
a numerical integration of a test satellite. Fig. 9 shows a clear de- 
pendence of e max with inclination. In a later section we will study 
and discuss the time evolution of these initial conditions to unravel 
their chaotic and quasi-periodic nature in detail presenting a proof - 
of-concept of detecting chaotic initial conditions. In this work we 
will discuss the location of mean motion resonances and postpone 
the study and discussion of other features (boxes located at polar 
and prograde orbits) in a future work. 

6.2 Fine structure of mean-motion resonances 

Fig. 7b shows details of the retrograde region with a <G [0.14; 0.20] 
AU ([0.39,0.56] R H , [293,419] R Jup ) and/G [130°; 180°]. Several 
vertical structures are observed corresponding to the location of 
chaotic mean-motion resonances. If n and nj denotes the mean- 
motion of the satellite and Jupiter respectively, then the semi-major 
axis of a satellite in a n : nj mean-motion resonance with Jupiter is 
given by 



nj 



-2/3 / v -1/3 

/ mj + Mq \ 

\ mj J 



(14) 



where aj,mj denote Jupiter's semi-major axis and mass, respec- 
tively. Mq is the mass of the Sun. In the figure we show the location 
of several mean-motion resonances by arrows. See Table 1 for their 
nominal locations. Several high-order mean-motion resonances are 
detected. 

When compared to the current population of observed irreg- 
ular satellites it is interesting to note the large scatter in (a J) el- 
ements of the Pasiphae group. Fig. 7b strongly implies that sev- 
eral members of this group are strongly affected by the dynamics 
within high-order mean-motion resonances. The dynamical con- 
sequences of mean-motion resonances on the orbits of irregular 



satellites were already reported by Saha & Tremaine (1993) dis- 
cussing the n — 6nj ~ resonance of Sinope and possibly S/2001 
Jl 1 (Nesvorny et al. 2003). Showing a more compact distribution of 
the osculating elements members of the Carme family are on less 
inclined orbits / « 165°. It appears that the dynamical effects of 
the mean-motion resonances occurs at a smaller magnitude for this 
group possibly having the effect of a smaller dispersion in orbital 
elements. In addition the Ananke group shows also a small scat- 
ter in their orbital elements possibly due to the close proximity of 
the 7:1 mean-motion resonance. The coincidence between the scat- 
ter of orbital elements of retrograde satellites and the presence of 
mean-motion resonances is probably not by chance. We stress that 
this work does not address the dynamical significance and effects of 
mean-motion resonances. Our results raises the question of the dy- 
namical effects of mean-motion resonances on the orbital elements 
on a compact group of satellites. We will address this question in 
a future work and study the distribution of orbital elements of an 
initial compact group by gravitational scattering in mean-motion 
resonances. 



6.3 Chaotic regions and Kozai mechanism 

From the MEGNO map (Fig. 7a), we observe a clear distinc- 
tion between quasi-periodic and chaotic phase space regions. The 
general chaotic region for prograde satellites starts from a = 
0.14 AU (0.39 Rh,293 R Jup ) and outwards. Stable quasi-periodic 
orbits for the retrograde satellites are found for orbits with semi- 
major axis up to 0.195 AU (0.55 R H ,409 Rj up ). It is interest- 
ing to note that the chaotic phase- space of prograde satellites is 
larger when compared to the retrograde satellites. This asymme- 
try in (a, I) space was already discussed by Nesvorny et al. (2003) 
pointing out the difference between prograde and retrograde sta- 
bility limits. At larger semi-major axis the retrograde satellites are 
on dynamically more stable orbits when compared to the prograde 
irregular satellites. In this work the calculated MEGNO map con- 
firms this stability asymmetry which is mainly explained by the ex- 
istence of the evection resonance for prograde satellites (Nesvorny 
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Figure 10. Secular time evolution of osculating elements from integrat- 
ing IC-I. Elements a,e were smoothed using a window width of w x AT — 
2500 x 40 days = 274 years. For the inclination / we used w = 5000. The 
black vertical bars shows half the window width. The elements 03, X denotes 
the (retrograde) longitude of pericenter and mean logitude of the satellite, 
respectively. P measures the orbital period of the satellite and P is the 
orbital period of the Sun in the jovicentric system. (j>i = 5X Q — X — 403© . 
<|>2 = 5X & - X- 203© - 203. (|>3 = 5^0 - X — G5 — 05 — 2£l & . (j)4 = 5^ - X + 
203-6^0. 
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Figure 11. Polar representation of the temporary circulation and librations 
of 05 — 05© about 0° (05 — 05/ = 180°) in the 5:1 mean-motion resonance of 
IC-I. The outer radius corresponds to e w 0.33 and the inner radius corre- 
sponds to e « 0.41. The alternation between libration and circulation is a 
clear sign of chaotic dynamics. 



et al. 2003; Yokoyama et al. 2008, and references therein). In geo- 
metrical terms this resonance locks the satellites apocenter towards 
the direction of the Sun and maintains this alignment for a period 
of time. In this orbital configuration solar perturbations accumulate 
and the satellite eventually escapes from Jupiter when 05 — X® ~ 0. 
An interesting feature shown in the map is the chaotic 'horizon- 
tal cone' for polar orbits at around / ~ 90° extending from a = 
0.04 AU (0.11 R H ,84 R Jup ) to a = 0.12 AU (0.34 R H ,251 R Jup ). 
This region was already studied intensively by Carruba et al. 
(2002); Nesvorny et al. (2003) showing that test satellites with in- 
clinations in the range 70° < / < 110° are short-lived orbits with 
survival times less than 10 7 years. Analytical perturbation theory 
(Carruba et al. 2002) showed that the accumulation of secular solar 
perturbations is the driving force causing the excitation of 
orbit eccentricity to large values. This mechanism is known as the 
Kozai resonance. Depending on the initial value of 05 circular or- 
bits can reach high eccentricities and collisions with the inner reg- 
ular satellites (or with Jupiter itself) are expected to occur. Thus 
the Kozai resonance provides an efficient dynamical mechanism 
to remove an initial population of near-polar irregular satellites. In 
near-polar orbits when studying Fig. 9 we notice that the Kozai 
regime is well identified with large eccentricities attained in the 
range 40° < / < 140°. However, it needs to be mentioned that the 
eccentricity excitations might strongly depend on the initial 05 (Car- 
ruba et al. 2002). 



6.4 Chaoticity versus quasi-periodicity 

Previously we have reported on the presence of chaotic dynamics 
as detected from calculating 7, (7) on a grid of initial conditions 
in (a,/)-space. In the following we will study in detail two initial 
conditions that were detected to be either quasi-periodic or chaotic. 
Fig. 7 marks the two initial conditions by dots (surrounded by cir- 



12 Hinse et al. 



Initial conditions 


ao (AU) 


eo 


/o(°) 


Region 


IC-I 


0.173915 


0.20 


153.968 


CH 


IC-II 


0.176179 


0.20 


153.968 


QP 



Table 2. Initial conditions for the two test satellites IC-I and IC-II in and 
close to the 5:1 mean-motion resonance as shown in Fig. 7. The remaining 
Kepler elements were set to co = 0°, £2 = 0° and M = 90°. CH and QP 
denotes initial conditions in the chaotic region and quasi-periodic region, 
respectively. 



cles) with IC-I and IC-II. Both initial conditions have the same orbit 
inclinations with different semi-major axis. We refer to Table 2 for 
details on numerical values in the initial conditions. We have then 
numerically integrated both orbits using the Radau integrator with 
initial step size of 0.01 days and tolerance parameter of 10~ 13. Cal- 
culations were done in double precision enabling the 'high' output 
precision in MERCURY. Initial conditions for the Sun were ob- 
tained from the JPL Horizon Ephemeris. We sampled osculating 
elements every 40 days. To maintain consistency we integrated the 
orbits over 5 x 10 3 Pj years. For both orbits the maximum relative 
energy error at the end of integration were smaller than 10 -13 . 

As suggested from MEGNO initial condition IC-I exhibits 
chaotic behaviour. We present the results of our single orbit cal- 
culation in Fig. 10 and Fig. 11. We have obtained those figures by 
successively applying the running window time average smoothing 
technique as outlined previously. In addition to the orbital elements 
(a,e,I) we also plot the time variation of 05 — G5 along with four 
critical angles (|>i , (|>2 , <|>3 , <|>4 (see details in the figure caption). In the 
bottom panel of Fig. 10 we also show the time evolution of 7, (Y) 
as a function of time. We obtained this plot from a GBS integration. 
Quantitatively at time 60000 yrs (Y) « 18 and a visual inspection 
shows a clear trend of divergence of (Y) over time. The presence 
of chaos is qualitatively best seen in the time evolution of the angle 

05 — G5 . In Fig. 1 1 we give a polar representation of this angle. A 
similar approach was also adopted by Saha & Tremaine (1993). In 
the beginning this angle librates around 05 — G5 = 0°. After approx- 
imately 30000 years this libration mode switches into circulation 
for a short time period and then returns to the libration mode. At 
the end the angle circulates. This qualitative change between dif- 
ferent modes of motion is characteristic of motion near a separatrix 
and hence chaotic motion is concluded. In the polar representation 
times of temporary librations are shown as 'banana' -shape curves 
and circulations are indicated by full circles. 

A more elongated banana corresponds to a larger libration am- 
plitude of 05 — G5 . From Fig. 10 it is apparent that the libration am- 
plitude changes with time. The initial conditions are chosen with 
the orbital inclination to be initially outside the Kozai resonance 
for which co of the satellite starts to librate around either 90° or 
270° (Carruba et al. 2002; Nesvorny et al. 2003; Yokoyama et al. 
2003). No large eccentricity variations are expected. It is impor- 
tant to note the difference of the location of the libration centre 
found for the test satellite started at IC-I when compared to the 
libration behaviour of this angle for several major satellites. Saha 

6 Tremaine (1993); Whipple & Shelus (1993) report that for the 
retrograde satellites Pasiphae and Sinope the angle 05 — 05/ librates 
about 180°. The difference from this work is the choice in reference 
system. In a heliocentric system the angle 05 — G5 for IC-I would 
librate about —180° which is consistent with the general trend as 
reported in Saha & Tremaine (1993). 

Since IC-I is near the 5:1 mean-motion resonance (compare 
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Figure 12. Time evolution of a test satellite with initial condition IC-II. The 
same window width were used as in Fig. 10. Note: The range in semi-major 
axis, eccentricity and inclination is smaller than in Fig. 10. This time the 
angle 05 — G5 is circulating as opposed to the chaotic case. 
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Figure 13. Polar representation of the {a, I)— space of Jovian irregular satellites for two different initial orbit eccentricities. Left panel: eo = 0.10. Right panel: 
eo — 0.20. In both panels the MEGNO indicator (Y) is color coded from 1.5 to 4.0 (see the electronic version for colors). The initial semi-major axis (ao) is 
given in units of Jupiters Hill radius Rh, and Io denotes the initial orbit inclination. The diagonal lines represents orbit inclinations at 45° (prograde satellites) 
and 135° (retrograde satellites). The central vertical line represents the polar orbit at Iq = 90°. For comparison black dots indicate the initial conditions studied 
in (Carruba et al. 2002). 



with Table 1) we have systematically explored and looked for the 
existence of critical angles of the form Murray & Dermott (2001); 
Morbidelli (2002); Yokoyama et al. (2003) 



(15) 



with Y,ki i = and the sum of the coefficients of the nodes being even 
as is required by the d'Alembert's rules. The strength of a given 
critical angle depends on the power in eccentricity. From Fig. 10 
we have a max then (|>i , (|>2, <|>3 , <h circulates prograde sense. 

Also whenever I is minimum (maximum) then e is maximum 
(minimum). Also whenever librates around 71 then 7 is at a min- 
imum and constant (at time 35000 years). We also see that at the 
end of the time evolution of the satellite the angles §i,§2 and §3 
are switching between libration around 71 and prograde circulation. 
At some times we see a temporary resonance lock in §2 m the 5:1 
mean-motion resonance. 

Turning our attention to the time evolution of the quasi- 
periodic orbit started at IC-II we report on the following results. 
Fig. 12 shows the smoothed time evolution of the osculating ele- 
ments (a,e,I) along with 05 — 05 © . In the bottom panel plot we again 
plot the time evolution of 7, (Y) over 60000 years. The smoothing 
follows the running window time average techniques as outlined 
previously. Contrary to the chaotic orbit the angle 05 — G5 now cir- 
ulates only and the secular time variation of the semi-major axis, 
eccentricity and inclination are characterised by quasi-periodic os- 
cillations. 

Based on the preceeding comparative study we conclude that 
MEGNO is a reliable numerical tool for detecting chaotic dynam- 
ics on short time scales. We plan to use this tool in future work 
on the dynamics of irregular satellites for the major planets in the 
Solar System. A particular interesting subject of study would be 
the change of mass of Jupiter and the corresponding change in the 
topology structure of phase- space at a given time during the growth 
phase of Jupiter. 



6.5 Comparing with previous work 

In the following we compared our MEGNO maps with numerical 
studies published previously in the literature. Fig. 13 shows a polar 
representation of Fig. 7a for two different values of initial eccentric- 
ities eo = 0.10 and eo = 0.20. In Fig. 13a,b the final value of (Y) 
after 60000 yrs of integration time is color coded with yellow indi- 
cating chaotic dynamical time evolution and blue indicating initial 
conditions exhibiting quasi-periodic dynamics. The dynamical and 
collisional evolution of irregular satellites and their lifetimes has 
been studied by Carruba et al. (2002) and Nesvorny et al. (2003). 
We have compared Fig. 13 with the results of long-term integra- 
tions of hypothetical irregular satellites conducted by Carruba et al. 
(2002, fig. 8, fig. 9). Similar studies can be found in (Nesvorny et al. 
2003, fig. 9). Initial conditions of test satellites studied in Carruba 
et al. (2002) have been superimposed by black dots in both maps 
following the exact array of initial conditions as presented in Car- 
ruba et al. (2002). The initial semi-major axis ranges from 0.08 AU 
to 0.20 AU with spacing Aa = 0.02 AU. The initial inclination is 
35° to 70° for the prograde and 110° to 145° for retrograde satel- 
lites with spacing Ai = 5°. In the following discussion it is impor- 
tant to stress the difference in the dynamical models used. In this 
work we only consider the Sun- Jupiter-test particle system without 
considering collisions with other bodies or ejections. Carruba et al. 
(2002) includes the perturbative effects of the major planets and 
includes collision and ejection criteria within Jupiter's Hill sphere. 

When comparing our results for prograde test satellites (eo = 
0.10) with (Carruba et al. 2002, fig. 8) all test particles with life- 
times less than 10 Myrs (ao ^ 0.14 AU) are located in (or close 
to the onset of) the chaotic region shown in Fig. 13a. Initial con- 
ditions started in the quasi-periodic region have stable orbits over 
1 Gyrs for low-inclination (7 < 50°, a < 0.12 AU) orbits to 10 
Myrs for high-inclination orbits (Iq ^ 65°, ao ^ 0.08 AU). It is im- 
portant to note that although our MEGNO map shown in Fig. 13a 
indicates quasi-periodic dynamics for the high-inclination orbits 
(for ao = 0.08 AU to 0.12 AU) those initial conditions are short- 
lived due to the presence of the Kozai cycle which opens up a 
route for those satellites to reach into the region of the orbits of 
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the Gallilean satellites. As demonstrated by Carruba et al. (2002) 
such high-inclination test satellites will experience close encoun- 
ters or collisions with the massive regular satellites of Jupiter. For 
the retrograde orbits (with eo = 0.10) a comparison allows the con- 
clusion that all test satellites started in (or close to) the chaotic 
region have lifetimes smaller than 1 Gyrs and all test satellites 
with initial conditions in the quasi-periodic region have lifetimes 
over 1 Gyr with the exception of high-inclination orbits exhibiting 
Kozai cycles in their orbital eccentricities and inclinations. Simi- 
lar conclusion are obtained when comparing the results shown in 
Fig. 13b (e = 0.20) with (Carruba et al. 2002, fig. 8). A final inter- 
esting point to mention is the qualitative difference in (a, I) -phase 
space topology of irregular satellites when varying the initial ec- 
centricity. The prograde quasi-periodic region seen in Fig. 13a at 
/o ^ 35° with ^o/Rh x cos(/o) ~ 0.4 has significantly decreased in 
Fig. 13b when eo = 0.20. Furthermore the 5:1 MMR for retrograde 
orbits manifests itself more prominently in Fig. 13b when com- 
pared to Fig. 13a. This suggests that the structure of (a,/)-phase 
space topology is strongly dependent on initial eccentricity and a 
survey of this region of parameter space is currently ongoing. 



7 DISCUSSION AND CONCLUSIONS 

We have introduced, described and applied the MEGNO chaos in- 
dicator to the dynamics of jovian irregular satellites. Our results are 
based on the elliptic restricted three-body problem considering a 
test particle (irregular satellite) in a jovian planetocentric reference 
system perturbed by the Sun on an elliptic orbit. Initial conditions 
for the numerical integrations of the orbits have been obtained from 
the JPL Horizon Ephemeris Generator. Basic concepts and proper- 
ties of MEGNO are reviewed and described as well as details on the 
practical computation of Y(t) and (Y)(t) based on the variational 
equations. Numerical tests (see Fig. 2) have been carried out to de- 
tect and avoid artificial numerical chaos arising from the natural 
time discretisation of the applied numerical integration algorithms. 
Our test determined an optimal choice in the absolute and relative 
error tolerances required in the Gragg-Bulirsh-Stoer algorithm to 
detect real chaotic dynamics. We have calculated the MEGNO fac- 
tor for several known irregular satellites for the prograde and ret- 
rograde cases. Our calculation suggests that prograde satellites are 
more chaotic in contrast to retrograde orbits (see Fig. 4) though 
the detected chaos is very weak. To clarify this suggestion more 
detailed computations including also additional planetary pertur- 
bations along with Solar tides are necessary. Such calculations re- 
sulting in an estimate of the Lyapunov indicator are currently in 
progress for the major irregular satellites. 

It is important to note that every numerical tool capable of 
distinguishing between quasi-periodic and chaotic dynamics has 
limitations with regards to claiming quasi-periodicity considering 
only a limited period of time. This certainly is also the case for 
the MEGNO technique. In Fig. 3 we computed Y(t) and (Y)(t) 
for the retrograde irregular satellite Ananke. The plot shows that 
after 1 Myr the orbit of Ananke exhibits a chaotic orbit with 
(F)(l Myr) « 4.0. In the 60 kyr integrations the orbit of Ananke 
would possibly be interpreted as quasi-periodic with (Y)(t) deviat- 
ing only slightly from 2.0. 

Considering 35100 orbits we calculated the MEGNO indica- 
tor on a large grid in (a, I) -space known to be occupied by observed 
irregular satellites (Fig. 7). The resulting map revealed several in- 
teresting dynamical structures and we compared our results with 
previous studies addressing the question of the orbital stability of 



jovian test satellites. We found good qualitative agreement between 
chaotic (quasi-periodic) and unstable (stable) regions as was found 
previously (Carruba et al. 2002; Nesvorny et al. 2003). In particu- 
lar we confirm the asymmetry of the stable region when compar- 
ing prograde and retrograde satellite orbits Nesvorny et al. (2003). 
Retrograde satellite orbits have access to a larger volume of phase 
space characterised by orbit stability. This result is in contrast to the 
prograde satellites for which chaotic orbits and its associated insta- 
bility occurs (for e = 0.20) at a « 0.13 AU (0.37 R H ,272 R Jup ) 
and onwards (cf. Fig. 7). A small region of (a,/)-space at a w 0.14 
AU have been detected to indicate quasi-periodicity. In addition 
we detected the presence of mean motion resonances of retrograde 
satellites with the Sun. The location of several high order mean 
motion resonances were determined and compared with the present 
population of retrograde irregular satellite families. We find that the 
orbital elements of the members of the Pasiphae family are largely 
scattered as opposed to the Carme group exhibiting to occupy a 
more compact region in (a,/)-space. We preliminary explain this 
excess in scatter of the osculating elements due to the close prox- 
imity to mean-motion resonances. This postulate will be subject to 
a separate study currently ongoing addressing the question whether 
high-order (retrograde) mean motion resonances are capable of dis- 
persing an initial compact group of satellite members. 

To support our results obtained from MEGNO we also cal- 
culated and compared two initial conditions associated to two ret- 
rograde satellite orbits. The first were chosen to be close to the 
5:1 mean-motion resonance and the second initial condition were 
chosen to be just outside the location of this resonance. We then 
searched for signs of chaoticity and quasi-periodicity to validate 
the results obtained from calculating MEGNO. Applying a time- 
running smoothing window on the osculating elements our anal- 
ysis of the single orbit computations support the results obtained 
from MEGNO. Initial conditions started in chaotic regions are as- 
sociated to libration/circulation of resonant angles. Quasi-periodic 
initial conditions show only circulating behaviour. 

Motivated by the success of applying the MEGNO technique 
to the dynamics of irregular satellites we plan to conduct a large 
parameter survey to identify further chaotic regions within the Hill 
sphere of Jupiter. In addition we plan to include giant planet per- 
turbations and generate similar MEGNO maps of observed pop- 
ulations of irregular satellites in orbit around the remaining giant 
planets in the Solar System. 
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